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We present a set of benchmark calculations for the Kohn-Sham elastic transmission function 
of five representative single-molecule junctions. The transmission functions are calculated using 
two different density functional theory (DFT) methods, namely an ultrasoft pseudopotential plane 
wave code in combination with maximally localized Wannier functions, and the norm-conserving 
pseudopotential code Siesta which applies an atomic orbital basis set. For all systems we find 
that the Siesta transmission functions converge toward the plane-wave result as the Siesta basis is 
enlarged. Overall, we find that an atomic basis with double-zeta and polarization is sufficient (and 
in some cases even necessary) to ensure quantitative agreement with the plane-wave calculation. 
We observe a systematic down shift of the Siesta transmission functions relative to the plane-wave 
results. The effect diminishes as the atomic orbital basis is enlarged, however, the convergence can 
be rather slow. 

PACS numbers: 72.10-d,73.40.Cg,73.63.Rt 



I. INTRODUCTION 

First-principles calculations of electrical conductance 
in nano-scale contacts represents a main challenge in 
computational nanophysics. The interest for this type 
of calculations began in the mid-nineties where advances 
in experimental techniques made it possible to contact 
individual molecules thereby making it possible study 
the transport of electrons through true nano-scale struc- 
tures [HQ. Apart from the scientific interest, the devel- 
opment of reliable simulation tools for nano-scale quan- 
tum transport is relevant in relation to the continued 
miniaturization of conventional semi-conductor electron- 
ics, but also for the introduction of the new generation 
of molecule based electronics. 

It has by now become standard to calculate conduc- 
tance in nano-scale contacts by employing a combina- 
tion of non-equilibrium Green's function theory (NEGF) 
and ground state density functional theory (DFT). The 
resulting NEGF-DFT formalism provides a numerically 
efficient way of evaluating the Landauer-Biittiker con- 
ductance due to electrons moving in the effective Kohn- 
Sham (KS) potential without having to calculate the 
scattering states explicitly. It has been applied exten- 
sively to a number of different systems ranging from pure 
metallic contacts, over organic molecules to carbon nan- 
otubes suspended between metallic electrodes. Overall 
the approach has been successful in describing quali- 
tative features and trends [U, 0], however, quantitative 
agreement with experiments has mainly been obtained 
for strongly coupled systems such as metallic point con- 
tacts, monatomic chains, as well asiunctions containing 
small chemisorbed molecules @, H, 0| ■ 

It is generally accepted that the NEGF-DFT method 
only provides an approximation to the true conductance 
- even if the exact exchange-correlation (xc-) functional 
could be used, and the quality of the result is expected 
to be strongly system dependent. Moreover, it is not easy 



to estimate the effect of using approximate xc-functionals 
such as the LDA or GGA. We mention here that more so- 
phisticated methods for quantum transport based on con- 
figuration interaction, the GW method, time-dependent 
DFT, and the Kubo formula have recently been pro- 
posed [1, H, EU EH, E2- However, these schemes are 
considerably more demanding than the NEGF-DFT and 
at present they cannot replace NEGF-DFT in practical 
applications. 

Irrespective of the validity of the NEGF-DFT approach 
and the role played by the approximate functionals, it 
remains important to establish a general consensus con- 
cerning the exact result of a NEGF-DFT calculation for 
a given xc-functional and specified system geometry, i.e. 
a benchmark. Although this might seem trivial, the 
present situation is rather unsatisfactory as different re- 
sults have been published by different groups for the same 
or very similar systems (several examples will be given 
in the text). Perhaps the best example is provided by 
benzene di-thiolate between gold contacts where the cal- 
culated conductance vary with u p to a two orders of mag- 
nitude for similar geometries @, lH, EI EE El, El El • 

The relatively large variation of the results indicates 
that the conductance, or more generally the elastic trans- 
mission function, is a highly sensitive quantity. Indeed, 
the implementation of the open boundary conditions 
defining the scattering problem represents some numer- 
ical challenges. Indeed, small errors in the description 
of the coupling between the finite scattering region and 
the infinite leads as well as improper fc-point samplings 
in supercell approaches, can introduce significant errors 
in the resulting transmission function. 

In this paper we take a first step towards establish- 
ing a common reference for NEGF-DFT calculations by 
performing benchmark calculations for a set of five rep- 
resentative nano-scale contacts. The benchmarking is 
achieved by comparing the transmission function ob- 
tained using two different and independent, albeit simi- 
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lar, NEGF-DFT methods: In one case the Hamiltonian 
is obtained from the Siesta DFT program which uses 
a basis of localized pseudo atomic orbitals (PAOs) to- 
gether with norm conserving pseudopotentials. The sec- 
ond method applies a basis of maximally localized Wan- 
nier functions (WFs) obtained from the Dacapo DFT 
code which uses plane- waves and ultra soft pseudopoten- 
tials. In both cases we use periodic boundary conditions 
in the directions perpendicular to the transport direction 
and we apply the PBE xc- functional [l9[ . 

The five reference systems we have chosen for our 
benchmark study are: (i) A monatomic gold chain with 
a single CO molecule adsorbed, (ii) A 3-atom Pt chain 
suspended between Pt electrodes, (iii) An H 2 molecule 
bridging two Pt electrodes. (iv) Benzene-dithiolatc 
(BDT) between Au electrodes, and (v) Bipyridine be- 
tween Au electrodes. The systems have been chosen ac- 
cording to the criterion that both experimental data as 
well as previous NEGF-DFT calculations are available 
in the literature. Furthermore they are representative 
in the sense that they cover a broad class of systems: 
homogeneous and heterogeneous, computationally sim- 
ple (one-dimensional) and more complex, and strongly 
as well as weakly coupled. 




FIG. 1: (color online). Deviation between the WF and Siesta 
transmission functions for for the five reference systems stud- 
ied. The dashed line indicate zero deviation from the WF 
transmission. Notice that the Siesta results converge toward 
the WF result as the PAO basis is enlarged. 

A main results of our work is summarized in Fig. [T] 
where we show the overall deviation 
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between the transmission functions calculated using the 
WF and PAO basis sets, respectively. The cutoff energy 
E is taken to be the energy above which the WFs are 
no longer able to reproduce the exact KS eigenstates of 
the system which is typically ~ 3 eV above the Fermi 



level. For all the systems we find that the deviation A 
decreases as the Siesta basis is enlarged meaning that the 
Siesta transmission functions converge toward the WF 
result. We take this as evidence for the correctness of 
the WF results and the justification for the use of the 
term benchmark calculation. 

In general we find that the double-zeta polarized 
(DZP) basis provides very good agreement with the WF 
basis, whereas the single-zeta polarized (SZP) and, in 
particular, the single-zeta (SZ) basis can produce sub- 
stantially incorrect features in the transmission function. 

The paper is organized as follows. In Sec. [Til we briefly 
review the NEGF-DFT formalism and introduce the two 
specific implementations used in the present study. In 
Sees. IIIUIVII we present the benchmark calculations for 
the five reference systems and in Sec. I Villi we give our 
conclusions. 



II. METHOD 

In this section we first outline the NEGF-DFT method 
which has become standard for nanoscale conductance 
calculations. The two specific NEGF-DFT implemen- 
tations applied in the present work are then introduced 
and their key parameters are discussed. We then consider 
the important issue, common to both methods, of how to 
treat periodic boundary conditions in the plane perpen- 
dicular to the transport direction. We end the section 
with a discussion of the advantages and disadvantages of 
the two methods. 



A. NEGF-DFT 

The zero temperature, linear response conductance 
due to non-interacting electrons scattering off a central 
region (C) connected to thermal reservoirs via two bal- 
listic leads (L and R), can be written as 



G = GqT(sf), 



(2) 



where T(e) is the elastic transmission function and Go = 
2e 2 /h is the quantum unit of conductance. Using the 
NEGF formalism, Meir and Wingreen have derived a 
very useful formula which expresses the transmission 
function in terms of the Green's function of the central 
region H, 



T(e)=Tr[G r (e)T L (e)G a (e)T R 



(£)]• 



(3) 



In this expression the trace runs over the central region 
basis functions and Tl/r. is obtained from the lead self- 
energies (defined in Eq. §5§ below) as T L / R = i(Y< L / R — 

In the NEGF-DFT method both the leads and cen- 
tral region are modeled by the effective KS Hamilto- 
nian, /jks = _ 5 V 2 + v c ff(r). The self-consistent effec- 
tive potential consists of the well known parts v c g = 
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v C xt + vh + v xc . Introducing a basis of localized orbitals, 
{<pi}, we define the Hamiltonian and overlap matrices by 

Hij = {^i\h K s\^j) and = {<j>i\4>j), respectively. In 
the original derivation by Meir and Wingreen the basis 
was assumed orthogonal, but the generalization to non- 
orthogonal basis sets shows that Eq. © still holds when 
the Green's function is defined as [21| 

G{z) = [xS c -H c -^l{z)-^r{z)]-' l . (4) 

Here the matrices He and Sc are the blocks of H and S 
corresponding to the central region basis functions. The 
retarded Green's function, G r (e), is obtained for z = 
e+iO + , and the advanced Green's function for z = e — iO + 
or G a = (G r )t. 

The self-energy of lead a is given by 

T, a (z) = {zS Ca - H Ca )g a {z){zS aC - H aC ), (5) 

where Hc a and Sc a are the coupling- and overlap matrix 
between basis functions in the central region and lead a, 
respectively. g° is the surface Green's function describing 
the isolated semi-infinite lead, g^(z) = [zS a — Ha] -1 , 
which can be calculated recursively using the decimation 
technique [22|. 

B. Method 1: Wannier functions from plane- wave 

DFT 

In method 1 the Kohn-Sham Hamiltonian is ob- 
tained from an accurate plane-wave pseudopotential 
DFT code [23j . The ion cores are replaced by ultra- 
soft pseudopotentials [H| and we use an energy cutoff 
of 25 Ry for the plane wave expansion. The Kohn-Sham 
eigenstates are transformed into partly occupied Wannier 
functions (WFs) [2|| which are used to obtain a tight- 
binding like representation of the Hamiltonian. The WFs 
are constructed such that any eigenstate below a selected 
energy, Eq, can be exactly represented by a linear com- 
bination of WFs. In the applications we have chosen Eq 
in the range of 2-4 eV above the Fermi level. In this way 
the accuracy of the plane- wave calculation is carried over 
to the WF basis for all energies relevant for transport. 

By performing separate DFT calculations for the (peri- 
odic) leads and C we obtain a set of WFs for each region. 
Note, that C always contains a few buffer layers of the 
lead material on both sides of the nano-contact to en- 
sure that the KS potential at the end planes of C has 
converged to its value in the leads. Since the WFs in 
the lead in general will differ from those in the outer- 
most lead unit cells of the central region, care must be 
taken to evaluate the coupling and overlap matrices Hc a 
and Ssa- Notice also that although the WFs by construc- 
tion are orthogonal within each region, WFs belonging to 
different regions will in general be non-orthogonal. For 
more details on the construction of the WFs and the 
calculations of the Hamiltonian matrix for the combined 
L — C — R system we refer to Ref. [26|. We shall refer to 
the results obtained from method 1 as the WF results. 



C. Method 2: PAO Siesta basis 

Method 2 is based on the DFT code Siesta H3 which 
uses finite range pseudoatomic orbitals (PAO) [28| as 
basis functions and Troullier-Martins norm conserving 
pseudopotentials [2{J. As in method 1, the Hamiltonians 
for the leads and the central region are obtained from sep- 
arate calculations. Because the KS potential to the left 
and right of C , by definition has converged to the value 
in the leads, we can take the coupling between central 
region and lead a, Hc a , from the pure lead calculation. 
Note that this is in contrast to method 1, where the dif- 
ferent shape of the WFs in the periodic lead and the lead 
part of the central region makes it essential to evaluate 
the coupling matrix directly. Note also that this approx- 
imation, i.e. the use of the intra-lead coupling matrix 
elements (H aa ) in Hc a , can be controlled by including a 
larger portion of the lead in C . In practice we find that 
3-4 atomic layers must be included in C on both sides of 
the junction to obtain converged conductances. 

In the present study we restrict ourself to the standard 
PAOs for Siesta: SZ, SZP and DZP. For the confinement 
energy, determining the range of PAOs, we use 0.01 Ry 
and for the meshcutoff we use 200 Ry. 

D. Common ingredients 

In both methods 1 and 2 we treat exchange and corre- 
lation effects with the PBE energy functional 19]. Fur- 
thermore, we impose periodic boundary conditions in the 
surface plane directions. This means that we are in fact 
considering the conductance of a periodic array of junc- 
tions instead of just a single junction. Instead of the 
localized basis function (f> n (r) (this could be a WF or a 
PAO) we thus consider the Bloch function 

X»k„ =]>> lk 'r R ^„(r-R||), (6) 
R 

where Rm runs over supercells in the surface plane and 
k|| is a wave- vector in the corresponding two-dimensional 
Brillouin zone (BZ). Since k|| is a good quantum num- 
ber, we can construct the Hamiltonian, i?(kii), for each 
k-point separately. This in turn implies that the conduc- 
tance per junction is given by the average 

G = 5>(k||)G(k||), (7) 

k|| 

where u>(k||) are symmetry determined weight factors. 
Unless stated otherwise, we have used a 4 x 4 Monkhorst- 
Pack k||-point sampling of the surface BZ, which for all 
the systems studied yields conductances converged to 
within a few percent [2(| [3(| • We take the Fermi level of 
the bulk lead as the common Fermi level of the combined 
system by shifting the levels in the central region by a 
constant. This is done by adding to He the matrix SS C , 
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where S = [Hl]o,o — [Hc]o,o an d the (0, 0) element cor- 
responds to the onsite energy of a basis function located 
near the interface between L and C. 

The main advantages of method 1 are: (i) The ac- 
curacy of the plane wave calculation carries over to the 
WF basis set. (ii) The WFs basis set is truly minimal 
and often results in even fewer basis functions than a 
SZ basis. The WF basis thus combines high accuracy 
with high efficiency. The price one has to pay is that the 
actual construction of well localized WFs is not always 
straightforward, and requires some user interaction - in 
particular for metallic systems. Also the lack of finite 
support of the WFs is unwanted in the context of trans- 
port, although in practice it is not a serious problem since 
the WFs are well localized. Finally, as already explained 
above, the risk of obtaining different WFs for two similar 
but non-identical systems renders it less straightforward 
to patch the parts together using the Hamiltonians ob- 
tained for the separate calculations. 

Most of the disadvantages of the WF basis are re- 
solved by the PAO basis set: by construction they have 
finite support and are identical as long as the atomic 
species on which they are located are the same. This 
renders it straightforward to patch together Hamiltoni- 
ans for separate sub systems as long as the KS poten- 
tial can be smoothly matched at the interfaces. On the 
other hand, to obtain an accuracy matching the WFs, 
one needs to use a significantly larger number of orbitals 
and thus longer computation times as compared to the 
WF method. 



III. GOLD CHAIN WITH CO 

In this section we calculate the conductance of an in- 
finite gold chain with a single CO molecule adsorbed. 
Scanning tunneling microscope (STM) experiments sug- 
gest that CO strongly depresses the transport of elec- 
trons through gold wires [3 11 ] - This has been supported 
by NEGF-DFT calculations [HJ which shows that the 
transmission function indeed drops to zero at the Fermi 
level. The use of infinite gold chains as leads is clearly 
an oversimplification of the real situation, however, the 
model seems to capture the essential physics, i.e. the 
suppression of the conductance, and furthermore is well 
suited as a benchmark system due to its simplicity. 

The geometry of the system is shown in Fig. ^a.). We 
use a supercell with transverse dimensions 12A x 12A 
and take all bond lengths from Ref . l32l: c?Au-Au = 2.9 A, 
d\u-c = 1-96 A, and dc-o = 1-15 A. The Au atom hold- 
ing the CO is shifted towards CO by 0.2 A. In method 1 
we obtain six WFs per Au atom and seven WFs for the 
CO molecular states. Due to the elongated bond length 
of the Au-wire, we found it necessary in method 2 to in- 
crease the range of the Au PAOs in order to converge the 
band-structure of the Au-wire. The confinement energy 
was therefore in this case set to 10 _4 Ry. 

In Fig. [DJb) we show the calculated transmission func- 




FIG. 2: (color online), (a) Central region used to describe a 
single CO molecule adsorbed on a monatomic Au wire, (b) 
Transmission functions for the Au wire CO system calculated 
using method 1 (WF) and method 2 for three different PAO 
basis sets. The transmission function at the Fermi level is 
indicated in the parenthesis following the legends. 



tion using three different PAO basis sets and the WF 
basis set. Overall, the PAO result approaches the WF 
result as the basis set is enlarged. For the largest PAO 
basis (DZP) the agreement is in fact very satisfactory 
given the differences in the underlying DFT methods, 
e.g. ultrasoft- versus norm-conserving pseudo potentials. 
The remaining difference can be further reduced by a 
rigid shift of the DZP transmission by about 0.15 eV. 

All transmission functions feature an anti-resonance 
near the Fermi level. However, upon enlarging the PAO 
basis the position of the anti-resonance shifts as fol- 
lows: (SZ) -0.27 eV, (SZP) -0.16 eV, (DZP) -0.06 eV, 
and (WF) 0.12 eV. Note that the position of the anti- 
resonance obtained with the WFs is approached as the 
PAO basis set is increased. Also, the curvature of the 
anti-resonance is improved as the PAO basis set is en- 
larged. The improvement in these features are, however, 
not directly reflected in the conductances indicated in the 
parenthesis following the legends in Fig. [2jb). The rea- 
son for the this apparent disagreement is the rigid shift 
between the PAO and WF transmission functions. 

We observe that our results differs from the calculation 
in Ref. [32|: While the latter finds two peaks in the energy 
range — 2 eV our converged transmission function shows 
a single broad peak. In general, both our PAO and WF 
based transmission functions present less structure than 
the transmission function reported in Ref.[32|. We suspect 
that these differences are related to the way the coupling 
V a c is calculated. 
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IV. PT CONTACT 

Atomic point contacts formed from late transition met- 
als such as Au, Pt, and Pd show very stable and repro- 
ducible features in conductance measurements (jj. This, 
together with the simplicity implied by their homogene- 
ity, makes them ideal as benchmark systems for trans- 
port calculations. Here we consider the conductance of 
a pure Pt contact for which both experimental conduc- 
tance data [33l l34l . |35| as well as theoretical calcula- 
tions [1, 0, l3d| are available. 

Conductance histograms obtained from mechanically 
controlled break junction experiments on pure Pt sam- 
ples show a peak near 1.5 Go, indicating that as a Pt 
contact is pulled structures with a conductance around 
1.5 G are frequently formed. NEGF-DFT calculations 
have shown that (zig-zag) monatomic Pt chains indeed 
have a conductance close to this value 0, 0, Hfl . More- 
over, the calculations predict an increasing conductance 
as the Pt chain is stretched and evolves from a zig-zag to 
a linear configuration. This effect has also been observed 
experimentally [35| . 




e-e F (eV) 

FIG. 3: (color online), (a) Supercell used for the DFT cal- 
culation of a short linear Pt chain between Pt(lll) surfaces, 
(b) The calculated transmission function using method 1 and 
method 2. The transmission at the Fermi level is indicated 
in the parenthesis following the legends. In the inset we show 
the transmission function in the important region near the 
Fermi level for the DZP basis set and the WF basis set. 

In Fig. [3Ka) we show the supercell used to model the 
scattering region of the Pt contact. The Pt contact is 
modeled by two four-atom pyramids attached to (111) 
surfaces containing 3x3 atoms in the surface plane. In 
order to ensure that the effective KS potential has con- 
verged to its bulk value at the end planes of the supercell 
we include 3-4 atomic layers (ABC-CABC stacking) on 
either side of the pyramids. The chain is formed by in- 
serting a single Pt atom between the apex atoms of the 



two pyramids. We have relaxed the contact region (pyra- 
mids+chain) while keeping the rest of the structure fixed 
in the bulk configuration with lattice constant 3.93 A 
and a distance of 14.60 A between the (111) surfaces. 
The cutoff energy used in the construction of WFs was 
set to £i? + 4.0 eV ensuring that the KS eigenstates below 
this value are exactly reproducible in terms of the WFs. 

In Fig. 0(b) we show the calculated transmission func- 
tions using method 1 and method 2. The qualitative 
agreement between the two methods is striking, however, 
only the SZP and DZP basis sets provide quantitative 
agreement with the WF result. The SZ basis set results 
in a down shift of the peak at —6 eV and an incorrect de- 
scription of the features in the important region near the 
Fermi level. Here the converged transmission function 
shows two peaks positioned at ep — 0.8 eV and just be- 
low the Fermi level, respectively. The main peak astride 
the Fermi level in fact consists of three smaller peaks, 
as seen more clearly in the inset for the DZP and WF 
basis set. These particular features in the transmission 
function were also observed in Ref. for a similar Pt 
contact, employing a method based on quantum chem- 
istry software and a description of the bulk electrodes by 
a semi empirical tight-binding Hamiltonian on a Bethe 
lattice [33| ■ Also, the calculated conductance of 2.3 Go is 
in agreement with our results, considering the structural 
differences. 

In Fig. [4] we show the calculated conductance of the 
Pt contact for three electrode displacements: 13.62 A, 
14.60 A, and 14.75 A. The three configurations corre- 
spond to an unstrained Pt chain, the chain just before it 
snaps, and the broken chain, respectively. 
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FIG. 4: (color online). Conductance for three different con- 
figurations during the stretching of a small Pt chain. Config- 
urations 1,2, and 3 correspond to the unstrained chain, max- 
imally strained chain, and a broken chain, respectively. The 
contact atoms are shown in the insets 

All basis sets, except for the SZ, are able to reproduce 
the trend of increasing conductance prior to rupture. The 
SZ basis set underestimates the absolute conductance by 
nearly 0.5 Go in the strained and broken configurations as 
compared to the WF result. The conductance calculated 
with the SZP and DZP basis set is almost identical and 
shows results more consistent with the WF basis for all 
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three configurations. 



V. PT-H 2 -PT CONTACT 

In this section we consider the simplest possible molec- 
ular junctions, namely a single hydrogen molecule be- 
tween metallic Pt contacts. Like the metallic point con- 
tacts, the Pt-H2-Pt junction shows stable and repro- 
ducible behavior in conductance measurements. In par- 
ticular, a very pronounced peak close to IGo appears 
in the conductance histogram obtained when a Pt con- 
tact is broken in a hydrogen atmosphere [33]. Although 
reported conductance calculations show significant vari- 
ation (see below), there have been given substantial ev- 
idence that the structure responsible for the peak con- 
sists of a single hydrogen molecule bridging the Pt con- 
tacts m, ni]. 

Several groups have published NEGF-DFT calcula- 
tions for the transmission function of the Pt-H2-Pt sys- 
tem. Most find a conductance of (0.9 — 1.0)Go, but also 
much lower values of (0.2 — 0.5)Go have been reported. 




FIG. 5: (color online), (a) Supercell used to model the central 
region of the Pt — H2 — Pt junction, (b) Transmission function 
for the Pt Hydrogen bridge. The transmission function at 
the Fermi level is indicated in the parenthesis following the 
legends. 

For the benchmark calculations we use the same setup 
as in Sec. IIVI with the central Pt atom replaced by a 
hydrogen molecule, see Fig. EJa). The relevant bond 
length determining the structure after relaxation of the 
Pt pyramids and the hydrogen atoms are dp t _H = 1.7 A 
and g?h-h = 1.0 A. 

In Fig. [5Jb) we show the calculated transmission func- 
tions. Like in the case of the Pt contact, the agreement 
between the different calculations is striking, especially 



in the important region around the Fermi level. The SZ 
basis set reproduces the qualitative features of the larger 
basis sets, but introduces a considerable down shift of the 
low-lying peaks. 

The very good agreement between the two methods 
indicates that the transmission function for this system 
is rather insensitive to the basis set. We stress, however, 
that a proper k|| -point sampling of the transmission func- 
tion is crucial to obtain meaningful results independently 
of the quality of the basis set. Restricting the calculation 
to the r point yields a transmission function with a (un- 
physical) peak at the Fermi level [g]. We note in passing 
that such a peak is present in the transmission function 
reported in Ref. 0. Such unphysical features resulting 
from an insufficient k|| -point sampling are not proper- 
ties of the molecular junction, but are rather due to van 
Hove singularities in the quasi one-dimensional leads [3(| • 
The results reported in Ref. E3 are based on Siesta DFT 
code and show good agreement with our results. The 
conductance obtained in one of the early theoretical cal- 
culations [391 ] on the hydrogen molecular bridge are con- 
siderably lower than our and most other results. The 
calculational method applied in Ref. |H is, however, the 
same as applied in the study of pure Pt contacts [36| . 
which agrees well with our results as discussed in SeclIVI 
We speculate if this could be related to the smaller size 
of the Pt cluster used to model the electrodes in [3^ as 
compared to[36|. Another possibility for the discrepancies 
is the use of the B3LYP energy functional in [3^ instead 
of an LDA/GGA functional used in most other works. 



VI. BENZENE-l,4-DITHIOL (BDT) BETWEEN 
AU(lll) SURFACES 

The Benzene- 1,4-dithiol (BDT) molecule suspended 
between gold electrodes was among the first single- 
molecule junctions to be studied and has become the 
standard reference for calculations of nano-scale conduc- 
tance. Depending on the experimental setup, measured 
conductances vary between 10~ 4 Go and 10 _1 Go jH, 
EH El EH El] > while the calculated value s typica lly lie in 
the range (0.05 - 0.4) G @, El M El EES El El- 
In general it has been found that the transmission func- 
tion is strongly dependent on the bonding site of the S 
atom [l8|, l47| . while variations in the Au-S bond length 
only affects the transmission function weakly [46|. 

As our objective is to establish a computational bench- 
mark for the Au-BDT system we choose the simple junc- 
tion geometry shown in Fig. [Ha). The S atoms are 
placed in fee hollow sites of the Au(lll) surface and the 
molecule has been relaxed while keeping the Au atoms 
fixed in the bulk crystal structure. We use an Au lat- 
tice constant of 4.18 A, and a distance between the 
Au(lll) surfaces of 9.68 A. With these constrains the 
relevant bond lengths are: g?au-s=2-45 A, g?s-c = 1-73 A, 
and e?c-H = 1.09 A. 

In Fig. [3b) we show the calculated transmission func- 
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FIG. 6: (color online), (a) Supercell used to model the central 
region of the Au(lll)-BDT-Au(lll) system with S at the fee 
hollow site, (b) The calculated transmission functions. Note, 
that the SZ transmission function has been omitted for clarity. 
The transmission function at the Fermi level is indicated in 
the parenthesis following the legends. 

tions (the SZ result has been omitted for clarity). Notice 
that we plot the transmission function only up to 2 eV 
above the Fermi level. This is because the the WF trans- 
mission at larger energies is sensitive to the parameters 
used in the construction of the WFs, in particular the 
cutoff energy Eq, and thus we cannot be sure about the 
WF result above 2 cV + Ep. 

The three transmission functions agree very well in the 
energy range from 2 eV below the Fermi level to 1 eV 
above the Fermi level, while only the DZP result agrees 
quantitatively with the WF result in the entire energy 
range. We again notice the down shift of the PAO trans- 
mission functions relative to the WF result. 

The presence of a broad transmission peak positioned 
~ 1 eV below the Fermi level is in qualitative agreement 
with previous results @, [IS EE Ea El, For more 
stretched configurations, i.e. for larger values of the S-C 
bond length, than the one used in the present study, the 
broad peak splits into two more narrow peaks J26| . 

The transmission function presented in Ref. y was ob- 
tained using a method very similar to our method 2, how- 
ever, the reported conductance of 0.4 Go is almost twice 
as high as our DZP results of 0.24 Go- The large conduc- 
tance arises because the transmission peak closest to the 
Fermi level is considerably broader than what we find. 
If, however, we restrict the ki| to the T-point we find the 
same broadening as in Ref. y and a very similar conduc- 
tance of 0.37 Go- Another feature of the r-point only 
transmission function is that the second peak positioned 
at ~ 3 eV below the Fermi level separates into a number 
of more narrow peaks. 



In Ref. Il6| the transmission function is calculated from 
the LMTO-ASA method and averaged over 36 irreducible 
k||-points. Both the width and the position of the two 
peaks in the transmission function at 1 eV and 3 eV below 
the Fermi level, are in good agreement with our results. 
The height of the former peak is, however, lower than 
in our calculation and this reduces the conductance to a 
value of 0.07 Go- We suspect that this difference could 
be due to differences in the adopted contact geometries. 



VII. BIPYRIDINE BETWEEN AU(lll) 
SURFACES 

As the last reference system we consider a bipyridine 
molecule attached between two gold-electrodes. STM 
experiments on bipyridine molecules in a toluene solu- 
tion [5l| show that the conductance of Au-bipyridine 
junctions is quantized in multiples of O.OlGo which was 
interpreted as the formation of stable contacts containing 
one or more molecules. The conductance is expected to 
be sensitive to the details of the contact geometry [52j . 
however, for the benchmark calculation we choose the 
simple case of binding at the on-top site of a fiat Au(lll) 
surface, as shown in Fig. 0[a). The Au electrodes are 
the same as used for the BDT molecule in Sec lVIl the 
Au(lll)-N distance is 2.180 A, while the electrode dis- 
placement is 11.53 A. 

The transmission functions calculated using either 
PAOs or WFs are shown in Fig. [3b). At first it is noted 
that the overall structures of the transmission functions 
are very similar. In the Siesta calculations, the position 
of the narrow LUMO peak which governs the transport is 
underestimated but converges towards the WF result as 
the PAO basis set is enlarged, see the inset of Fig. \T[b). 
The alignment of the LUMO energy level with respect 
to the Fermi level and its relation to charge transfer was 
studied in Ref. l53l . 

Several groups have investigated the transport prop- 
erties of bipyridine-gold junctions, and there is general 
agreement that the low bias conductance depends cru- 
cially on the details of the contact geometry. As different 
groups have chosen different geometries and models for 
the gold electrodes a direct comparison of the reported 
transmission functions is difficult. 

To the best of our knowledge the first theoretical paper 
on the bipyridine system is by Tada et al. (54|. In their 
calculations bipyridine is adsorbed on-top of an Au-atom 
of a rather small Au cluster, and the coupling to the infi- 
nite electrodes is modeled by a broadening parameter fit- 
ted to experimental data. The zero-voltage transmission 
function contains some of the same peak structures as 
we observe. Hou et al. [HI have published several papers 
on the gold-bipyridine juncition. Like Tada et al. they 
include only a few gold atoms in the ab-initio calcula- 
tion and treat the coupling to electrodes through a model 
self-energy term. The peak structure of the transmission 
function is quite different from ours. This could be due 
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to the small size of the gold clusters used to mimick the 
electrodes. While most other groups observe tunneling 
through the LUMO tail HH, [56] , Hou et al. argue that 
the transport is mainly taking place via the HOMO-2 
state. Calculations by Wu et al. [I3,[58| obtained using a 
Siesta-based transport code for bipyridine attached 
to the on-top site of a gold surface show overall goo d 
agreement with our results (see Fig. 7(a) in paper |57T). 
The minor differences are probably related to the fact 
that only the T-point has been used in the latter paper. 
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FIG. 7: (color online), (a) Supercell used to describe the 
central region of the bipyridine- Au(lll) junction, (b) Calcu- 
lated transmission functions (the SZ result has been omitted 
for clarity). The inset shows the dependence of the LUMO 
position on the basis sets. The transmission function at the 
Fermi level is indicated in the parenthesis following the leg- 
ends. 



VIII. CONCLUSIONS 



We have established a set of benchmark calculations 
for the Kohn-Sham(PBE) elastic transmission function 
of five representative single-molecule junctions using two 
different methods based on independent DFT codes: (i) 
A plane wave DFT code in combination with maximally 
localized Wannier functions, (ii) The Siesta program 
which applies finite range pseudoatomic orbitals. 

For all five systems we find that the Siesta result con- 
verges towards the WF result as the Siesta basis is en- 
larged from SZ to DZP with the latter yielding very good 
quantitative agreement with the WF transmission. In the 
Siesta calculations the transmission peaks relative to the 
peaks obtained with the plane-wave calculation are sys- 
tematically shifted toward lower energies. The problem 
can be overcome by enlarging the Siesta basis, however, 
the convergence can be rather slow. 
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